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Abstract. In computational relativity, critical behaviour near the black hole 
threshold has been studied numerically for several models in the last decade. In this 
paper we present a spatial Galerkin method, suitable for finding numerical solutions 
of the EiNSTEiN-DiRAC equations in spherically symmetric spacetime (in polar/areal 
coordinates). The method features exact conservation of the total electric charge and 
allows for a spatial mesh adaption based on physical arclength. Numerical experiments 
confirm excellent robustness and convergence properties of our approach. Hence, this 
new algorithm is well suited for studying critical behaviour. 
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1. Introduction 



1.1. Preface 

In 1998-1999 F. Finster, J. Smoller and S.T. Yau published static soliton- 
like solutions of the massive Einstein-Dirac equations in a spherically symmetric 
spacetime in ^ and D.W. Schaefer, D.A. Steck and J.F. Ventrella studied 
the critical collapse of that system numerically. Although the applied iterative Crank- 
NlCHOLSON finite differencing scheme does neither conserve the total electric charge 
exactly nor make use of adaptive mesh grid refinement methods, the authors were able 
to investigate critical collapse for low (m = 0.25) and moderate (m = 0.5) particle 
masses for suitable initial data (See jTTj for details). 

Yet, these non adaptive finite difference methods fail for simulations of critical collapse 
for large particle masses [m > 1) or propagating DiRAC waves as initial data. Therefore 
we performed a Galerkin discretization of the same equations and combined this with 
a physically motivated adaptive mesh grid refinement algorithm. In this paper we give 
a proof of the discretized conservation law, show some convergence tests and compare 
our results to those in [TTj. 

Finally we note, that critical behaviour of the massless Einstein-Dirac system 
has already been studied using iterative Crank-Nicholson finite differencing with 
adaptive mesh grid refinement by J.F. Ventrella and M.W. Choptuik in [T2j. 
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1.2. The basic model 



Throughout this paper we use Planckian units, where the speed of hght, Planck's 
action and Newton's gravitational constant satisfy c = h = Gn = 1- We consider a 
spherically symmetric spacetime in polar/areal coordinates {t,r,6,(p) with metric 

g = A{t, r) dt (g) dt - B{t, r) dr (g) dr - d6 ® d9 - sm^{0) dcf) ® d0, (1) 

where A{t,r) and B{t,r) are positive valued functions. We couple this metric by 
Einstein's equations to a massive 2-fermion system of two spin-^ particles in singlet 
state. In order to find a coordinate representation of the Einstein-Dirac equations, 
we constructed a complex 4x4 matrix representation of the Clifford bundle and an 
appropriate spin bundle with fibers isomorphic to C^, by using tetrad- formalism^ 
1^. This yields 7-matrices with {7'^,7'^} = 2 g^'^ 1, a representation of the (2,2)-signed 
semi-hermit ian scalar product (., .) on spinors and the connection coefficients of the 
Levi-Civita connection on the spin bundle. In analogy to the ansatz for the DiRAC 
spinors in flat spacetime interacting with a spherically symmetric potential, we write 
the wavefunctions in the form§ 



^1 
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a 



1 




1 




1 



a 




ism{9) sin(</))/5 
i cos(^) — sin(6') cos(0) ] /3 



i2a) 
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icos{9) + sin(6') cos(0)]/3 
— i sin(^) sin(0)/? 



where a{t,r) and (3{t,r) are complex valued functions. These wavefunctions satisfy 





47rr^v 













[i«r-i/?r] 



(3 a) 
(36) 



On each spacelike slice of constant t, denoted by S(t), we construct the fermion state 
of the two particles, given by 

1 



1^) 



V2 



1^1) ® 1^2) - 1^2) ® 1^1)] . 



(4) 



I We used Pauli-Dirac representation throughout this paper, but one ends up with exactly the same 
expressions for all relevant physical quantities by using any other representation within the same unitary 
class, e.g. Weyl's or Majorana's representations [TU] . 

§ The prefactors — r are slightly different from those used in |4ji to remove any time derivatives 

of A and B from DiRAC's equation and to make the normalization condition for the functions a and P 
as simple as possible. 
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Assuming each particle carries || an electric charge q, the components of the charge 
current density vector J = g X]p=i (^p' l^i'p) ~ J^^ti would be 



B 









(5) 



and the total electric charge on each slice sums up to 



Q:= / (J,N)dS = g / ^t^^ + 
which implies the normalization [4j condition^ 

[\a\^ + \P\^] dr=l. 



dJ: = 2q I [\a\^ + \P\^]dr = 2q, (6) 



(7) 



With the ansatz ()2a|26|l the equations of motion for the wavefunctions ijji and 11)2 are 
independent of each other, i.e. they satisfy DiRAC's equation separately: 

[iV -ml]'iljp = (8) 

Here m > is the particle mass and V := denotes the DiRAC operator 

constructed from the metric (jT]). The stress energy tensor of the system is just the 
sum of the stress energies of both wavefunctions, its components are given by 

. 2 



p=i 



From the ansatz ()2 al2 6|) we get 

1 [A 



Too 
Toi 

T22 



27rr2 V B 
1 



1 

27rr2 
1 



Im aa* + pp' 



Re a*P - a(3* 



(10) 



A 

— Im [aa* + = Tio 



27rr/B 



Re (a*/?' -a'/?*) 
Re {a* (3) and Tsa 



1 



27rr/B 



Re sin2(^). 



Plugging the ansatz ()2al26|) into (jH)) yields two first order evolution equations, which 
are linear in a and (3: 



a 



(3 ^/A[3 - im^fAa 

r 



'11a) 



II We don't couple the particles to the electromagnetic field in this paper. The particle charge q is 
introduced just for convenience. 

^ The same result can also be optained by normalizing the quantum states on each slice S(i), i.e. by 
requiring {ipp\'4>p) = (^|^) = 1- 



13= + I 



a 



a 



Aa + imVAp 



[lib) 



Einstein's equations involving the stress energy Q leads to the two ordinary 
differential equations in the radial variable 
1 

A ~ r 
1 

'B " r 



[B - 1 + 2[a* - a' f3* + a(3*' - a*' p]] 
[ 1 - S + 2 [ - a'P* + a[3*' - a*'/3 ] ] 



;i2a) 
(126) 



Arn , A , 



which could be identified as the Hamiltonian constraint and the slicing condition in 
spherical symmetry, and an evolution equation for B: 

l = -J^\a'a*-aa*' + l3'(3*-(3(3*'] (13) 
B r \ B 

The equations ()126|1 and ()13|1 are consistent. One easily checks that the evolution of 
B w.r.t. ()13|) yields the same as solving ()126|) on each slice S(t). So one of them is 
redundant. 



1.3. The evolution boundary value problem 

We are looking for regular and asymptotically flat solutions of ()llalll6|) and 
(jl2Jl2/)l13|l . Let T > 0, / := [0,r] and J := [0,oo). We seek solutions of (|1 1 all 1 6jl 
and p2all26ll3p . such that for all t G / the variables a and f3 as functions of r are in 
the SoBOLEV space iy^'°°(J) fl H^{J) and A and B as functions of r are elements of 
W^'°°{J). This leads us to the "boundary conditions" 

lim A{t, r) = lim B{t, r) = lim B{t, r) = l\fte I. (14) 

r— »oo r— +00 r^O 

Because of the factor ^ in the ansatz ()2 al2 6|) . a and (] have to fulfill 

a{t,r) = 0(r) and /3(t,r) = 0(r) Vt G / for r 0. (15) 

For the numerical treatment of our problem, and in order to guarantee the asymptotical 
flatness of the solution, we consider the DiRAC matter to be concentrated inside a ball of 
coordinate radius R > around the center of spherical symmetry. Due to Birkhoff's 
theorem, the part of spacetime at regions where r > R has to be Schwarzschild 
|13j . Therefore we consider the following boundary /initial value problem on J x J/j with 
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Jr := [0,R] for a, /3, A and B: 

PDEs: (lllalll/)!) 
ODEs: ^i'2aiV2t^ 

a{t, 0) = P{t, 0) = a{t, R) = P{t, R) 

B{t,0) 
A{t, R) 
a{0,r) 




1 



B{t,R) 

ao{r) 
Mr) 



for (t, r) G (0, T) x (0, R) 
for (t, r) el X (0, i?) 
for t G / 

for t G / (16) 
for t G / 
for r G Jr 
for r E Jr 



Here ao,/3o e n i/o(JK) denotes some suitable initial data. Because of the 

boundary conditions for a and (3, outgoing DiRAC waves are reflected at r = /?, which 
is an unphysical requirement. For initial data with a fast decay (i.e. exponential decay) 
for large r and by specifying R large enough, these reflected waves do not interfere with 
what is going on in the region of interest, that is, for small r. 



1.4- Conservation laws 

The total electric charge from ©, i.e. the normalization condition (|7j), and the total 
ADM-mass of spacetime are constants of motion. For the charge this can be seen by a 
straightforward computation jB] using ()llalll6|) : 

J poo POO 

d f . . r r. . ...1 . ^^^^ 



Q = J^I [\af + \P\ ] dr = 2Re 



aa* + (3(3* 



dr 



2 Re 



— I 



2Rez 



2Rez 



Re? 



2Rez 



/A i 
- a (3* + - 
t) z 



a(3* 



[a(3* + Q*(3] + zmVI [ \(3f - \af] 



dr 



[a(3* - (3a* 



dr 



^ [ a' (3* - P'a* + a(3*' - (3' a* - (3a*' ] 



[a'(3* - (3'a* -a(3*' + (3a*'] dr 



dr 



Re{a'(3* - (3'a*)dr = 



The conservation of the ADM-mass of spacetime is obvious from the following formula, 
which could be transformed into a conserved integral by using Einstein's equations 
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and P|: 



Madm(S) = - lim r 

/ r— >oo 



1 " 










^0 



Im ad* + 



2. Conservative discretization 

To find a conservative Galerkin discretization, we first consider the model problem 
a = - UP' - - I P - iQa tip (19) 



- tf(3'-'-fP~iga--hp 
2 r 



/3 = + if a' H — fa + ig(3 ha , 

2 r 

where f,g,h G C^{I,H^{Jr)) are real valued. We seek complex valued solutions 
a{t),l3{t) G Hq{Jr). In the sequel (., .) denotes the scalar product in L?{Jr). The 
weak formulation of ()19p reads 



2 r 



+ z/q;' H — /'a + z^r/? ha, w ] , 

2 r 



(20) 



for all v,w & Hl^Jpi). In the spirit of the method of lines, we first carry out the spatial 
Galerkin semi-discretization of ()20|) . It can be obtained by replacing the function 
space Hq{Jji) with a finite dimensional subspace Vn- Please note that these subspaces 
have to contain continuous functions that vanish for r = and r = R. As a basis 
for Vn we choose real valued global shape functions Cj G Hq^Jr) for i G {1, . . . ,n}, 
n := dimVn G N. Thus we can represent spatially discrete approximations of a and P 
in (1201) by 



«(f,r) ^ ^a'(t)ei(r) and P{t,r) ^ ^ P'{t) ei{r) . 



i=l 



i=l 



Next, we define the n x n GALERKlN-matrices E, F, G and H with components 



^ , := {gei,ej) and if^j := ^^/iei,ej^ 



From the integration by parts formula we get 



R 



dr = — 



R 



dr = -K 



and the above matrices obey the algebraic properties 

= E, = -F, = G and = H 



(21) 



(22) 



(23) 



(24) 



Then, we tackle discretization in time: We recall that (fTIH) has the character of a wave 
equation with a quadratic first integral, the total electric charge Q, see Sect. 11.41 Exact 
conservation of quadratic first integrals is guaranteed by the so-called implicit midpoint 
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rule, the simplest representative of Gauss collocation Runge-Kutta timestepping 
methods jHl Sect. IV.2.1]. Formally, this timestepping scheme is second order accurate 
with respect to the size of the timestep. To describe its application to the semi-discrete 
system, we fix two instances in time < < tm+i < T, introduce the (local) timestep 
r := tm+i — tm, the notations 

(25) 



Ctr. 



and analogous for (3m, A(3 and (3. The fully discrete version of then reads 
-EAa= - iPp- iGa - iH(3 

T 



(26) 



r 



EAp=+ iFa + iG(3 - iHa. 



Obviously, each timestep involves the solution of a linear system of equations. The 
discrete total electric charge is given by 

Qm = 2q[ alEam + PlEpm ] ■ (27) 

Because of the properties ()24|). it is conserved exactly along the time evolution generated 
by the discrete equations (jJEI): 
1 



2gr 



AQ = a'^AAa + P^AA(3 + i 2 Im {ai,Aam+i + PlA(3m+i) 



(28) 



:5GIR 



[ (3^ Fa - a^F(3 ] +i [ a'^H(3 + P^Ha ] +i [ /3^G(3 -a'^Ga] +iS 



eR 



eR 



eR 



0. 



We point out, that this conservation property holds for any choice of trial/test space 
Vn, as long as the same Vn is used for every timestep. 



3. Computer implementation 

3.1. Variables 

We implemented the discretization scheme described above for the six real variables 
Xa, Ya, Xb, Ife, a and b, defined by 

a = Xa + iYa, f3 = Xb + iYb, A = e'' andB = e\ (29) 

where the representation of A and B is motivated by the essential positivity of the 
metric coefficient and the presence of logarithmic derivatives in ()12aj) . ()12&|) and (fT^ . 
In these variables DiRAC's equations read 



a — b , X 

X, = + e— Y! + - 









e 2 













(30 a) 
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a — b 




"1 X. 


e 2 


Xt,-e2 




- 




r 


a — 5 






e 2 


V„ + et 








r 



a — tt 




e 2 









and Einstein's equations are 



and 



+ — [XI + Y^^ - XI - + ^et [X^X, + Y^Y, 



6 = ^ e^ [ X'Xa - XX + - XX : 



(30 
(30 c) 
(30 c?) 

(31a) 
(316) 



3.2. Initial conditions 

We implemented two classes of initial conditions. The one parameter family 
Gaussian's, which was also used in [TT], is defined by 



X,(0,r) : = 



S 2 re 



(32) 

(S) of 
(33) 



F,(0,r) :=X,(0,r) =n(0,r) = 0. 
Further we considered the 8 parameter family (S^, Sf,, Pa, Pb, "^a, ^6, ^«;,p): 

2nr^ 



Xa{0,r) 



Ya{0,r) 



n(o,r) 



1 p 

24~2 COS (\l/i 



- COS + 



24"2 sm (\1> 



sin \& 



r(p+i) E 

23-i sin (^^) 



- COS 



r sin *b + 



27rr 
27rr 
27rr 



P P 4e2 



(34) 



r^e ''^6 



r^e *^b. 



Since critical behaviour is always studied by using a one parameter family (see e.f 
one of the 8 parameters will be chosen to vary, while the other 7 remain fixed. 



0), 
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3. 3. Mesh and shape functions 

As Galerkin trial and test spaces for X^, Fq, X^., and we rely on spaces of 
piecewise polynomial continuous functions on a spatial mesh covering Jr. In other 
words, we use classical Lagrangian finite element spaces. More precisely, the mesh is 
a not necessarily equidistant grid on J^j, whose N cells, X G N, are the intervals 

Jn ■■= [r„_i,r„] with Tq = 0, = R, A„ := r„ - r„„i. (35) 

Following the finite element philosophy, we use bases of the piecewise polynomial 
function spaces that only comprise locally supported (global) basis functions. They 
are assembled from the following local shape functions (see figure [T]) : 

for r E Jn 



1 




Elir) 



otherwise 
(r - r„) + 1 

1 + 5nN 



(36) 







for r G (r„-i,r„ 

for r = Tn 
otherwise 



El{r) 



-r) + 1 

1 + ^nO 



Kir. 



Kir) 



4 
A 



2 ir - r„) {vr, 



64 

3A| 



for r e (r„_i,r„) 

for r = r„_i 
otherwise 

Xn - r) for r G J„ 
otherwise 

An 



(r - Tn) (r„ - A„ - r) r 







for r E Jn 



otherwise 



The parameters N, Tn and A„ characterizing the mesh can be redefined at every time 
step (See below for the details of our adaption algorithm). The global shape functions 
are obtained from the local shape functions through 

E TtE:,{r) (37) 



ef(r) 



!,n=l 



where T-^" are the combination coefficients 



rpOn 

i 

T, 



IS 

In 
- is 



6l6?+' + s'is: 



(3^ 



S I 
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Local shape functions: 




_-| 5 I I I I I I I I I I I 

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 

x:=(r-r 

Figure 1. 



Local shape functions of different polynomial degree on an interval of the spatial grid. 



S t 

-,B^ 



With this definitions the sets Sb '■= {ef } of global basis functions contain compactly 
supported functions of polynomial degree p = for B = and BN + 1 compactly 
supported functions of polynomial degrees p G {1, . . . , B} for B G {1, . . . , 3}. For the 
finite element space with basis S'b we adopt the notation Sb- This spaces will play the 
role of the abstract Galerkin trial/test space from sect.|21 

3.4- The discrete equations 

We represent the variables Xa, Ya, Xf, and 1^ of the DiRAC fields due to the discussion 
of section El as in ()2H1 by functions from Sb for some fixed parameter B G {1,...,3}. 
For the variables a and b we just use a linear approximation by functions in Si: 

N+l N+l 

a{tm, r)^Y^ oT^'e] and r) b'^'^e] (39) 

1=1 i=l 

The representations have to be compatible with the boundary conditions in (fTH|l . which 
implies 

-Ym,{BN+l) _ ^^m,{BN+l) _ ■5^m,(BAr+l) _ ^rm,{BN+l) _ ^ 

b""'^ = and a'"'^+^ = _5"^,jv+i^ 

We construct the discretization of a weak formulation of the equations ()30all32p 
analogous to ^I^ . whereby we have to use test functions from Sq for a stable 
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discretization of the ODEs ()Hla|H16j) . We define the coefficients 





:= (e' 


H ; 


- (e 




:= (el 


■Dijk 




Eij 


- ({ 



^ ^3 

Hijki 

Jijkl 



(IplpB B J}\ 
{ ,,.2 ^k-i^l ) 



(41) 



r 



and finally arrive at the following discretized versions of 

1 1 ^m,% rm,i m+l,j frn+l,i 

= -[ X-'^ - Xr'' ] A. + ^ e^ '—^^ [ + Y^'^^ ] S.,. 

7~ 



1 d"^' + a 



{ yr^' + yr-^'^' ] a,. + - e — ^ — [ y^'' + y^ ] d 



IJS 



o = ^[ Yr - y^'' e- — '—^^ — [ xr + xr''' ] B.,s 



^ — [^r-^' + xr''' ] c.,s - ^ e'^^ — [xr + ^r'" ] a,. 



= ^ [Xr- - Xr^'^ ] A. - ^ e^^ '—^^ [Yr + yr'^^ ] B.,s (42a) 



1 /I*"'* _|_ ^JTl+l,* 

[ y-'^- + ] c,,. + - e — [ Yr + yr''' ] n 



1 1 /I™'* ?)"^'* -J- rt^+li* trn,+l,i 

o = l[ Yr - yr'' ]As+l — '—^^ — [ xT^^ + xr'' ] 



= - a^'+^^'Eis + 
= - + 



+ 4m e 2 



1 n'"^''^ -4- 



Lm+l,i 

e — 1 



1 - e 



F,, + 4 + ] G,^.^ (425) 

F,, + 4 [Xr^'^xr^'-^' + ] G,,, (42c) 



Hijks 



+ 8e^~ 



= - [6^"'^ - 6"^+^'^] + e 



Hjks 



(42 d) 



+ Here we omitted the summation sign to shorten the notation. Analogous to Einstein's summation 
convention, we sum over the indices i, j and k where they appear twice. 
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Note that although we are free to decide which of the equations (jl2fejl or (fT!?|l we 
consider as redundant, the discrete systems ()42 al42 6142 cj) and ()42al426|42Q!j) are not 
exactly equivalent! 



3.5. The algorithms 

The following algorithms have been implemented in MATLAB: 

(i) An initial mesh can be given, initial conditions can be chosen from the families 

or (jH^ and the software projects them in L^(Jr) to a linear combination of 
the form ()21|) . All necessary integrals are evaluated by Gaussian quadrature [H], 
whose order can be specified in the parameter Fgo- 

(ii) The initial metric is obtained by solving the equations ()42 6p42 cjl implicitly by 
applying Newton's algorithm until the equations are satisfied up to a specified 
residuum PNtoi or a given maximum number FNmax of iterations is reached. 

(iii) Before each time step the mesh is updated due to the following procedure: The 
initial mesh structure is represented by a density function p(r), such that for all 
u,v e Jr 

/ pdr = Number of cells in ["Ujf ]. (43) 

J u 

The automatic mesh update ensures 



py B dr = Number of cells in [u,v] (44) 
by dividing a cell into two cells if 

I p\^dr^^ p>-^e^ + p^e^ > Path, (45) 

where Path is a threshold that can be specified as a parameter. The above criterion 
means the number of cells per physical arclength in radial direction is given by the 
initially specified function p on each spacelike slice S(t). This was motivated by 
the principle of equivalence. 

(iv) If the size r of the time steps is chosen to be too large, the discrete domain of 
dependence does not encompass that of the continuum system and the resulting 
time evolution may yield spurious solutions. Therefore we choose 

r := tjn+i -tm = Pxt ■ min l^Xje 2 j ^ (46) 

where P\r < 1 can be specified as a parameter. 

(v) For each time step the equations ()42ail424 are solved implicitly by applying 
Newton's algorithm until the residuum drops below a specified threshold PNtoi 
or the given maximum number PNmax of iterations is reached. In the latter case the 
time step is rejected if the residuum exceeds PNmaxtoi and repeated with a smaller 
Tnew := Par A ' ''old- After Pmaxtry refinements of T the program exits with an error 
message. 
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(vi) The residues of the equations ()42all42c|l and (j42djl as well as the change in the 
total electric charge and the ADM-mass are monitored and thus, can be used for 
accuracy checks. 



4. Numerical experiments 

4.1. Convergence Tests 

In order to study the convergence of the scheme, we conducted simulations on an 
equidistant static meshgrid, i.e. with automatic meshgrid and timestep refinement 
turned off. We performed 9 runs due to all possible combinations of the parameters 

S G {1, . . . , 3} and e {120, 240, 480} (47) 

for the initial data family ()33|) using the fixed parameter values 



-Pgo — 150, FNtol — 9-10 PNmaxtol — 10 PNmax — 30 

Path = 1.25, Pxr = 0.1, FarA = 0.5, F^axtry = 50, (48) 

m = 0.25, S = 0.3, = 5 and T = 3.125. 

Let M(_B,Ar)(r) be the numerical solution from the {B, A)-run for one of the variables in 
{Xa, Ya, Xh, Yh, a, b} at fixed time T and u{r) denotes the corresponding exact solution. 
Due to Richardson's extrapolation, we assume there is a function /^(r) and apB > 0, 
such that 

u = u^B,N) + Jb ■ [XnY'' , (49) 
where Xn '■= ^ is the grid spacing. For each run, we evaluated U(^b,n) on an 
equidistant test-grid with Nt ■= 10000 gridpoints, which yields the values |M(^''Ar)| 
for n G {1, . . . , Nt} and computed the approximation 



Pb ~ meanne{i,.. .,Nt} In 



(n.) (n) 

^(B,120) ~ '^(3,240) 



{n) (n) 

""(5,240) """(5,480) 



ln(2) . (50) 



The resulting values for pb are summarized in the following table: 



B Xa Ya Xb Yh a b mea.n{Xa,...,Yb} i^sanja^;,} mean all 

1 2.73 2.25 4.77 2.59 1.94 1.95 3.08 1.95 2.70 

2 1.95 2.13 1.96 1.90 3.13 3.05 1.98 3.09 2.35 

3 3.06 2.91 5.48 2.83 2.15 2.17 3.57 2.16 3.10 



4-2. Experiments near the blackhole treshold 

We tested our code for the initial data family ()33j) at particle mass m = 0.25 and 
compared our results to JT]. We specified the initial mesh by the cell distribution 
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function 

/ JcTT \ 

(51) 



D' := 1 + Pdc 



A;7r\ 
cos I -— 



from which the cell widths and the density function are computed by 

Afc := R ^'^ . and := ^ . (52) 

This results in an equidistant grid for Pqc = and for Pqc > the density of cells 
decreases with increasing r. For all runs we specified the following fixed parameter 
values: 

-Pgo = 150, PNtol = 9-10 PNmaxtol = 10 , -PNmax = 30 (53) 

-Path = 1-25, Pxr = 0.1, ParX = 0.5, Fmaxtry = 50, Pdc = 7. 

The remaining parameters to be specified for each run are m, R, B, Nq and S. 
For S small enough, the system collapsed to form a black hole, whereas all matter 
travels toward spatial infinity for large S. Using Crank-Nicholson finite differencing, 
SCHAEFER, Steck and Ventrella found the physical threshold Sth near Sthssv := 
0.412 for m = 0.25 in [H]. 

The decision, whether a black hole emerged, was based on the existence of a t and r, 
such that 

1 (54) 

r n{t,r) 

for a small parameter e^u > 0. By carefully varying S, we found the numerical blackhole 
threshold Sthz to satisfy 

0.41185 < Sthz < 0.41186 (55) 

for runs using B = 3 and A^o = 240. This confirms* the result of pi]. We show the six 
variables as a function of r at a fixed value of t at the end of the run with E = 0.41185 in 
figure El To investigate the details at the forming black hole horizon, we show a zoomed 
version in figure El In this run the black hole was detected at 

R 

Rbr ~ 0.052, Mbh := ~ 0.026 using esH := 0.9937. (56) 



The residue (difference of the right hand side from 0) of the redundant equation ()42dj) 
was saved for each global shape function at each timestep. The maximum over the 
global shape functions at each time is plotted in figure (H as well as the relative drift of 
the conserved total ADM-Mass, i.e. 

MADM(t) - Madm(O) 



ErrM(t) 



(57) 



Madm(O) 

* Since we don't know anything about the error with respect to the exact value of the treshold, 
we called Sthz the "numerical" treshold and publish all of its digits found. As mentioned above, 
ScHAEFER, Steck and Ventrella reported the critical value of Ethssv = 0.412, which was rounded 
to 3 significant digits, but was actually found to higher accuracy. Our result confirms at least this 3 
digits published in 
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Finally we show the development of the meshU due to our adaption algorithms in figure 
13 This consists of the number of timesteps performed to reach up to time t, the 
timestep width r(t) and the Number of cells N{t) as functions of t. By comparing these 
plots to their analogues for the run with E = 0.41186 in figures EHHl we conclude, that 
runs near the black hole threshold are well behaved and can be well studied before the 
numerics break down or are slowed down due to exorbitant mesh adaption. Therefore 
our algorithm seems to be suitable for studying critical collapse. 



tt Although we specified a value of A^o — 240, the number of cells in the plot starts at 299 because of 
the mesh adaption due to the initial metric, which takes place before the time development starts. 
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a(t,r)=XJt,r)+iYJt,r) and p(t,r)=X|^(t,r)+iY|^(t,r) at t=9.74 




A(t,r)=exp(a(t,r)) and B(t,r)=exp(b(t,r)) at t=9.74 



- A(0,r) 

- B(0,r) 

— A(t,r) 

— B(t,r) 



1 



10 



12 



Figure 2. 

Plot of the variables as functions of r for the run with i? = 12, i? = 3, A'o = 240 and E = 0.41185, 

after performing M = 5150 time steps. 
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Figure 3. 

Figure [3 zoomed to sliow tlie forming of the black hole. 
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Residues of the 3jb-equation: 

10- V;;.... . . . ^ ^ ^ ^ ., . . ^ ^ ^ ^ . . ., ^ ^ ^ ^ . . ^^^^ 




Change in the total ADM-Mass: 




Figure 4. 

Residues of the redundant equation 142 d|l and drift of the ADM-Mass for the run with i? = 12, _B = 3, 

iVn = 240 and S = 0.41185. 
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Time steps vs. time: m(t) 

6000 I 1 1 1 1 1 1 




01 23456789 

t 





Figure 5. 

Development of the mesh due to our adaption algorithms for the run with R = 12, B = 3, Nq = 240 

and S = 0.41185. 
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Figure 6. 

Plot of the variables as functions of r for the run with R — 12, B = 3, Nq — 240 and E = 0.41186, 

after performing M = 5950 time steps. 
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Residues of the 3 b-equation: 




Change in the total ADM-Mass: 

10"' I , 1 , 1 , 1 




10"^ I 1 1 1 1 1 1 1 1 1 1 

01 23456789 10 

t 



Figure 7. 

Residues of tlie redundant equation 142 d|l and drift of tlie ADM-Mass for the run with i? = 12, _B = 3, 

iVn = 240 and S = 0.41186. 
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Time steps vs. time: m(t) 




X 10 ' 



Time incrementations: T(t) 




Number of cells: N(t) 




Figure 8. 

Development of the mesh due to our adaption algorithms for the run with R = 12, B = 3, Nq = 240 

and S = 0.41186. 
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